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EXTREMA CLASSIFICATION 



Related Application 

[0001] This application claims the benefit of U.S. Provisional Application No. 60/461,782, 
5 filed April 10, 2003, incorporated herein by reference. 

Field of the Invention 

[0002] This invention relates to the field of seismic reservoir characterisation, and in 
particular to the automated mapping of structural geometries, like seismic interfaces and fault 
10 displacements. 

Background of the Invention 

[0003] Formation layers, lithological boundaries, sedimentary bedding, etc. in the 
underground can be defined through the interfaces between different lithologies, which 

15 produce seismic reflections due to impedance contrasts. These seismic reflections are referred 
to as seismic horizons, and interpretation of seismic horizons plays an important role in the 
structural characterization of 3D seismic data. Seismic horizons are commonly interpreted as 
being located along minimum, maximum, or zero crossing values in a seismic volume, and 
interpretations can be obtained by manual or automatic extraction of such surfaces. 

20 Structurally complex regions are a challenge to existing interpretation procedures, both 
automated and manual. For the manual interpreter, structural complexity produces 
ambiguities, and the interpretation process can become very time consuming. The continuity 
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constraint incorporated by most automated algorithms may fail to hold. In particular, 
automatic continuation of an interpretation across faults is a challenge. 

[0004] The present invention addresses this problem by laterally detecting and combining 
5 horizon segments having similar seismic waveforms. The mapping is not restricted to spatial 
continuity in the position of the seismic horizon, but instead determines possible lateral 
continuations based on similarity of signal shape. The inventive method represents 
waveforms as coefficients, typically seismic derivatives, that can be used to reconstruct the 
seismic trace in the vicinity of the extrema positions, such as by utilizing Taylor series 

10 expansions. The term "extrema" as used throughout this patent application includes any 
characteristics of seismic traces that can be used to track the positions of seismic horizons, 
such as minimum values, maximum values, zero-crossing values, midpoints between zero- 
crossing values and maximum or minimum values, etc. Derivatives that can be used in the 
Taylor series reconstruction and sub-sample accuracy extrema positions can be calculated 

15 based on orthogonal polynomial spectral decompositions of the seismic traces, such as 
described in U.S. Patent No. 6,240,370, issued May 29, 2001 to Lars Sonneland et al. 

[0005] The lateral continuations of horizons can be exploited to automatically estimate fault 
displacement, with a spatially varying displacement offset along the fault plane. Fault 
20 displacement influences the connectivity of hydrocarbon bearing lithologies, and assessment 
of the displacement provides an improved description of a reservoir. The present invention 
can calculate the displacement along pre-interpreted fault planes, for example fault planes 
that have been extracted from the seismic data using the methodology described in U.K. 
patent GB 2,375,448 B issued October 15, 2003 to Schlumberger Holdings Limited. 
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[0006] Accordingly, it is an object of the present invention to provide an improved method of 
processing and interpreting seismic data. In one embodiment, this involves extracting surface 
primitives, and grouping extrema positions where these surface primitives are both spatially 
continuous along the extrema of the seismic volume and are continuous in class index in the 
5 classification volume. 

Summary of the Invention 

[0007] Classification tools are widely used in both 3D and 4D reservoir characterization, for 
example in the mapping of 3D structures, lithological properties and production effects. In 
. 10 this invention, we extend the application area of the classification methodology into 
automated interpretation of seismic reflectors and fault displacement calculations. By 
classifying the seismic waveform along reflectors, we gain an improved automated 
interpretation that also performs well in structurally complex regions. The seismic waveform 
around the extrema positions can be represented by a set of coefficients that can be input into 
15 the classification process. These coefficients can also be used to reconstruct the seismic data 
waveform shape in the vicinity of the extrema positions, such as by using Taylor series 
reconstruction. One-point support for the reconstruction is an important element in the 
classification of seismic reflectors, as it allows the classification to be performed only along 
extrema positions while utilizing information regarding the waveforms in intervals around the 
20 extrema positions. This substantially reduces the number of data points to be classified, 
allowing the 3D classification to be run on a sparse 3D volume. The workflow of this 
preferred embodiment of the method, referred to as extrema classification, is outlined in 
Figure 1 . A related computer system and computer program product for implementing the 
inventive method are also described. 
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[0008] For existing fault surfaces, the present invention can automatically extract and match 
pairs of horizons on opposite sides of each fault. The displacement can then be calculated as 
the offset of these horizons in their intersections with the fault plane. This constitutes a fully 
automatic procedure for fault displacement assessment, providing displacement estimates that 
5 vary spatially along the fault plane. 

Brief Description of the Drawings 

[0009] Figure 1 is flow diagram outlining the extrema classification methodology. 

10 [0010] Figure 2 displays a vertical seismic section and the corresponding sparse extrema 
representation. 

[0011] Figure 3 illustrates the Taylor series reconstruction of a single seismic trace, 
displaying the observed trace, the reconstructed trace and the error. 

15 

[0012] Figure 4 shows the effect of increasing the number of terms in the Taylor series 
reconstruction of a seismic trace. 

[0013] Figure 5 illustrates how attribute vectors plotted in attribute space may be subdivided 
20 into different classes. 

[0014] Figure 6 shows how horizon patches may be extracted from an extrema classification 
cube, and how a horizon interpretation can be built. 
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[0015] Figure 7 shows how a extrema classification cube can be used to define seismic 
volumes. 

[0016] Figure 8 illustrates how the extrema classification technique can be used in the 
5 calculation of fault displacement estimates. 

[0017] Figure 9 exemplifies fault throw calculations along an extracted fault plane, using the 
extrema classification technique: fault surface, extracted horizon patches, and fault surface 
shaded by throw value. 

10 

[0018] Figure 10 gives an example of a horizon interpretation obtained in a faulted region. 

[0019] Figure 1 1 illustrates how extrema classification has been applied to extract a horizon 
in a region with weak seismic signal. 

15 

[0020] Figure 12 shows a fault population in which extrema classification has been used to 
estimate fault displacements. 

[0021] Figure 13 shows a schematic representation of a computer system and computer 
20 program product associated with the implementation of the inventive method. 

Detailed Description 

[0022] A flow diagram outlining the extrema classification methodology is shown in Figure 1 
as Extrema Classification Methodology 10. In this method, structural interpretations are 
25 placed on seismic extrema, and the first step in the workflow of the current invention is to 
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obtain a representation of extrema within the seismic volume, preferably with sub-sample 
precision. This is shown in Figure 1 as Generate Extrema Cubes 12 and involves identifying 
a plurality of extrema points associated with said seismic data, as discussed in more detail 
below. 

5 

[0023] Orthogonal polynomials may be used to reconstruct the seismic trace S(z) locally, 
such as by using the Volume Reflection Spectral decomposition technique (VRS) where 
orthogonal Chebyshev polynomials g, are used as the basis functions: 

S(z) = b 0 g 0 (z) + b x g x (z) + . . . + b n g„ (z) ( 1 ) 

10 This method is described in significantly more detail in U.S. Patent No. 6,240,370, issued 
May 29, 2001 to Lars S0nneland et al. 

[0024] Orthogonal polynomial reconstruction provides an analytic representation of the 
seismic signal from which high-accuracy positions of extrema can be calculated. The 

1 5 resulting extrema points may be represented through two sparse 3D volumes. A first cube can 
contain the amplitudes at the extrema, positioned at the voxels vertically closest to the 
extrema points. If the types of extrema being identified are zero-crossings, then a non-zero 
placeholder value (such as 1.0) could be used to mark the identified voxels. Voxels between 
extrema points are generally assigned zero value. A second cube contains values describing 

20 the vertical distances between the voxel positions where the amplitudes are stored, and the 
exact positions of the extrema, i.e., the sub-sample precision. The set of voxels containing 
extrema data is the same for the two cubes, but they contain amplitude and sub-sample 
position values respectively. Figure 2 shows a vertical seismic section 60 and a sparse 
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extrema representation 62 of this seismic section, a cross section of the cube containing these 
amplitude values. While the figures shown throughout this patent application are depicted in 
grey-scale, it will be understood that color displays of this information are customary and are 
preferable for many types of applications. 

5 

[0025] The second step in the workflow of the current invention is to generate a set of 
coefficients that characterize the seismic data in the vicinity of the extrema positions, which 
we refer to as "waveform attributes". The coefficients can be derivatives of the seismic 
signal at the single extrema data point and the seismic trace in the vicinity of the extrema 
10 point can be reconstructed using these coefficients in a Taylor series expansion: 

S(z + h) « S(z) + hS\z) + h 2 S 2 (z)/2\+... + h n S n (z)!n\ (2) 

This reconstruction gives a good fit to the seismic signal in a region around the point z, where 
the quality of the reconstruction depends on the number of coefficients/derivatives included 
in the Taylor series. Each trace in the 3D seismic volume is typically reconstructed using 

15 VRS decomposition, where the polynomial reconstruction (1) enables analytical calculations 
of the derivatives S*(z) at each sample point along the trace. Although multiple data points are 
involved in calculating the derivative values, the seismic signal in a region can afterwards be 
reproduced based on a set of derivatives obtained for one single data point. In particular, if 
derivatives are obtained only at extrema points along the seismic trace, this is sufficient to 

20 reconstruct the full seismic trace. This process is shown in Figure 1 as Generate Derivative 
Cubes 14. The derivatives are then extracted along the extrema positions and this is shown in 
Figure 1 as Extract Derivatives Along Extrema 16. Typically, a particular portion of the 
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subsurface area covered by the seismic data is selected by the user for more detailed analysis 
and this is shown in Figure 1 as Define Volume of Interest 18. 

[0026] Figure 3 shows an example of an observed seismic trace 70, a Taylor series 
5 reconstruction 72 of this seismic trace, and an error term 74 showing the difference between 
the observed seismic trace and the Taylor series reconstruction of the seismic trace. The error 
term is observed to be insignificant around the one-point support for the reconstruction. 
Figure 4 illustrates how the vertical region of insignificant error increases from left to right, 
as the number of derivatives included in the reconstruction (2) is increased. 

10 

[0027] Other types of reconstruction methods with one-point support can also be used in the 
inventive method. The one-point support of the reconstruction is; an important feature, since 
this enables the classification to be limited to extrema points only, but still be related to the 
overall shape of the seismic signal in the vicinity of the extrema points. 

15 

[0028] The next step involves forming groups of extrema points that have similar 
coefficients. Four alternative approaches to classification of seismic extrema may be utilized: 
three aiming at producing surface primitives for horizon interpretation or volume building; 
and one for fault displacement estimation. Common to all approaches is an underlying 
20 statistical model, assuming the derivative attribute vectors a k of length /?, where k is the voxel 
index, are Gaussian distributed within each class. The number of derivative attributes to 
include in the classification is typically selected by the user. Furthermore, the first derivative 
may be excluded from the classification, since they are zero, by definition, at maximum or 
minimum extrema points. 

25 
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[0029] The mean vectors // = (//, , . . . , fu n ) and covariance matrices £ = (£, , . . . , Z n ) of the 
Gaussian distributions describe respectively the position and shape in attribute space of 
classes 1 . The probability density function (pdf) of the complete set of K attributes 
vectors a = (a,,...,a^), given Gaussian parameters and class indexes c = (c, , . . . , c K ) of all 
voxels, is given as: 

/(a|c,//,S) = n/(a t \c k , Mck ,H Ck ), (3) 

Ar=l 



where the pdf of voxel k is given as: 

/(«, p^jjrop(-(«» "^)^ («, (4) 

[0030] The unknowns are the class indexes c k and the Gaussian parameters fi and S , which 
can be estimated using either supervised or unsupervised classification. The aim of the 
classification process is to subdivide the attribute vectors into n classes, providing a 
separation of the attribute space as illustrated for three attributes in Figure 5. Each shade in 
the right plot in the figure represents one class. The choice of an underlying Gaussian model 
can alternatively be generalized to other statistical models and the classification algorithm 
can be replaced by alternative classification approaches, such as neural networks. 

[0031] The volume of interest for performing the extrema classification can be chosen, for 
instance, as a vertical window with constant thickness or as the volume confined between two 
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pre-interpreted seismic horizons. This enables, for example, the extrema classification to be 
run only within a reservoir formation. 

[0032] In the first approach, supervised classification with picked seed points, the user picks 
5 seed points along reflectors of interest, each reflector representing a separate class. This is 
shown in Figure 1 as Pick Seed Points 20. Small, spatially continuous segments of the 
horizons are automatically extracted from the extrema cube, starting in the picked seed 
points, and the derivatives obtained along these extrema segments are applied as training 
data. The Gaussian attributes of each class are estimated from the training data, as the sample 
10 mean and sample covariance: 

where T c is the set of training data for class c and m c is the number of training data in T c . Two 
additional background classes are included: positive extrema not being training data; and 
1 5 negative extrema not being training data. Each extrema not being a training point is then 
assigned the class index maximizing the pdf (4) in that voxel. This process is shown in 
Figure 1 as Classify Supervised 22. 

[0033] In the second approach, unsupervised classification with a specified number of 
20 classes, the number of classes is provided as a user input (shown in Figure 1 as Select 

Number of Classes 24), and the attribute vectors are automatically divided into the specified 
number of classes, by estimating optimal class indexes and Gaussian parameters (shown in 
Figure 1 as Classify Unsupervised 26). The estimates are typically obtained by maximizing 
the likelihood function derived from the pdf (3). The maximum likelihood equations 
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constitute a non-linear system of equations, and cannot be solved analytically. Instead, an 
iterative scheme is applied to find the optimal parameter estimates. For an example of this 
type of interative scheme, see Celeux, G. and Govaert, G. (1995), Gaussian parsimonious 
clustering models, Pattern Recognition, 28(5), 781-793. 

5 

[0034] In the third approach, unsupervised classification with an unspecified number of 
classes, the algorithm starts with selecting at random a seed point among all extrema, and 
extracts from the extrema cube a small, spatially continuous horizon segment locally around 
the seed point. This process is shown in Figure 1 as Extract Training Horizon Segments 28. 

10 The extracted data constitutes a new class. The Gaussian parameters (5) are first estimated 
based on derivative attributes from this surface segment and for a background class of 
extrema not contained in the training segment. Next, the remaining extrema points are 
classified into either the class of the training segment or the background class, whichever 
maximizes the pdf (4). The initial horizon segment is then expanded along continuous voxels 

1 5 in the extrema cube, by including only those extrema points assigned to the same class. The 
procedure is repeated starting in new randomly selected extrema points not belonging to 
existing horizon segments, until all extrema points are contained in horizon segments with 
assigned classes. This produces a set of horizon segments, all with different class indexes. 
The number of classes is then reduced hierarchically, using classification to group horizon 

20 segments having similar seismic waveform shapes, i.e., similar average derivative attribute 
values along the surface segment. This process is shown in Figure 1 as Expand/Combine 
Using Classification 30. 

[0035] In a fourth approach, unsupervised classification for fault displacement estimation, 
25 fault surfaces are provided as input to the extrema classification algorithm. First, spatially 
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continuous horizon segments are automatically extracted from one side of the fault, each 
assigned to a separate class index. This is shown in Figure 1 as Extract Training Horizon 
Segments on One Side of Fault 44. Attribute values along these surface segments are then 
used as training data (5), and extrema points on the opposite side of the fault surface are 
5 assigned the class maximizing the pdf (4). Next, spatially continuous and class consistent 
horizon segments are automatically extracted from the classified extrema on the second side 
of the fault. The vertical sequence of the different classes on either side of the fault is 
compared, and classes being non-consistent within the vertical sequence across the fault are 
excluded from further calculations. This process is shown in Figure 1 as Classify on Other 
10 Side of Fault 46. 

[0036] As shown in Figure 1 , the inventive method may be used to interpret seismic 
horizons. The initial output from the extrema classification is a sparse 3D volume, with class 
indexes assigned to extrema points from the original input seismic volume. Thus, all 

15 classified data points are positioned on a reflector. Structural interpretations, e.g. horizon 
interpretations or volume interpretations, can be obtained by extracting and combining 
extrema points belonging to one, or a few, classes only. In order to obtain class consistent 
interpretations, small spatially continuous surface segments (patches) belonging to one class 
only can be extracted automatically from the extrema classification cube. This is shown in 

20 Figure 1 as Extract Horizon Segments 32. By filtering on properties like class index, 

position, attribute values, etc. attached to each patch, a set of patches can be combined into a 
final horizon interpretation. This process is shown in Figure 1 as Combine Segments Into 
Horizon Interpretation 34. In Figure 6, a sparse extrema classification cube is shaded by 
class index (upper left). Class consistent patches are extracted (upper right), and the number 
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of patches of interest reduced by filtering (lower left). The patches are combined into a 
horizon interpretation (lower right), shaded by its depth. 

[0037] A class consistent horizon interpretation, resulting from one class only, consists of 
5 extrema points with a similar shape of the seismic waveform in a neighborhood around the 
extrema, which are likely to belong to the same reflector. The building of the final horizon 
interpretations typically requires some manual effort. The extrema classification provides a 
set of patches, or geometry primitives, forming the basis for the interpretation. The 
classification index guides the interpreter in structurally complex regions, reducing the risk of 
10 interpreting along the wrong events. Overall, the extrema classification method can 

substantially improve the manual horizon interpretation process, easing the job and reducing 
the time spent by the interpreter. 

[0038] The inventive method may also be used to extract seismic volumes. When one or a 
15 few classes appear repeatedly vertically within a region, as illustrated in Figure 7, the seismic 
feature represented by these classes is better described as a volume than as a surface. A 
volume representation is obtained by defining the volumetric closure around the class or 
classes defining the seismic body, and should be interpreted as a seismic volume containing 
multiple reflectors of similar seismic response. This process is shown in Figure 1 as Extract 
20 Volumes 40. 

[0039] The inventive method may also be used in the estimation of fault displacements. 
Fault displacement is defined as the distance a horizon has been offset along the fault surface, 
due to faulting. Starting with an existing fault surface [see, for instance, Randen, T., 
25 Pedersen, S. L, and Sonneland, L. (2001), Automatic extraction of fault surfaces from three- 
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dimensional seismic data, In Expanded Abstr., Int. Mtg., Soc. Exploration Geophys., 551- 
554] (Figure 8a), the extrema classification provides a set of horizon segments located on 
either side of a fault surface. This is shown in Figure 1 as Extract Horizon Segments 48. 
Pairs of horizons belonging to the same class are identified on opposite sides of the fault, and 
5 the vertical sequence of classes is the same on either side of the fault (Figure 8b). A region 
around the fault surface, e.g., a fault zone, is excluded from the volume of interest, since the 
seismic signal may be corrupted by the fault in the close vicinity of the fault surface. 
Intersection lines between the horizons and the fault surface are obtained by extrapolating the 
horizon segments in their respective dip directions into the fault surface, represented by a 
1 0 plane placed optimally through the potentially non-planer surface (Figure 8c). 

[0040] Fault displacement calculations can be performed between the two intersection lines 
corresponding to a pair of horizons. The calculations are performed in the dip-slip direction, 
assuming no lateral movement along the fault surface. Displacements are calculated between 

1 5 pairs of intersection points, one from each of the two horizon segments, in a vertical plane 
placed perpendicular to the horizontal direction of the fault plane (Figure 8d). Each pair of 
horizon segments provides a displacement estimate varying laterally along a curve in the fault 
plane, and the full set of horizon segments provides a vertically varying set of estimates along 
such curves. A full coverage displacement estimate is obtained by interpolating between the 

20 local estimates along the curves. This procedure is shown in Figure 1 as Calculate Fault 
Displacement 50. 

[0041] The sample units of a seismic volume typically differs vertically and laterally, e.g., 
milliseconds in the vertical direction and meters in the horizontal direction. To avoid 
25 confusion between different measurement units, the displacement may be decomposed into 
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the vertical throw component and the horizontal heave component. Figure 9 illustrates a fault 
plane provided as input to the extrema classification procedure (top), the set of classified 
horizon segment pairs assigned different shades for each pair (middle), and the spatially 
varying throw estimate obtained from the horizon segments (bottom). 

5 

[0042] The primary advantage of the present invention over existing methodologies for 
automated seismic interpretation lies in its ability to track reflectors laterally in structurally 
complex regions, for example across faults. Figure 10 shows an example where the 
unsupervised extrema classification method has been applied to automatically extract a 
10 horizon that is cut by numerous faults. The reflector interpretation is shown both in map view 
and plotted on an inline and crossline section. The mapped reflector has a polygonal fault 
structure, and is interpreted automatically across a number of faults. Void regions in the 
interior of the fault blocks correspond to minor classification errors, and can easily be 
interpolated to produce a full coverage interpretation. 

15 

[0043] Figure 1 1 shows an example where the major challenge in horizon interpretation is 
the weak seismic signal. The reflector of interest is located in a relatively homogeneous sand 
formation, containing numerous thin, horizontal shale layers not thick enough to produce 
strong acoustic impedance contrasts. Extrema classification is applied to ease the horizon 
20 interpretation. The resulting interpretation has a good lateral coverage in most regions, but 
also shows some void regions. According to the extrema classification results, the void 
regions correspond to regions where the seismic response is lacking, or differs substantially 
from the rest of the reflector. This is valuable information, which can be related to structures 
like faults, channels, etc. 
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[0044] In Figure 12, the extrema classification technique has been applied to estimate fault 
displacements for a set of pre-interpreted faults. Figure 12 shows the throw component of the 
spatially varying displacement estimates. For many of the faults, the throw is observed to be 
at a maximum, at the center of the fault, decreasing towards the boundary, which is a typical 
5 geological feature of faults. 

[0045] The inventive method will typically be implemented on a computer system of the type 
shown in schematic format in Figure 13. The processes of the inventive method are 
translated into computer readable program code means and are recorded on a computer 

1 0 usable medium, such as Program Carrier 1 00, to produce a computer program product for 
processing and interpreting seismic data. This computer usable medium may consist of 
magnetic or optical recording media, such as floppy disks, hard drives, flash memory, CD- 
ROMs, or magnetic tapes. The contents of the computer program product are typically 
loaded onto a Storage Device 102 attached to a CPU 104 that executes the computer program 

1 5 recorded on the computer usable medium. The computer program may be loaded onto the 
Storage Device 102 either by placing the Program Carrier 100 into a reading device 
connected to the CPU 104 or by downloading the computer program using a network, such as 
a local area network or the Internet. When equipped with an operator input device (such as a 
Keyboard 106) and an output device (such as Display 108), the computer system shown in 

20 Figure 1 3 is capable of processing and interpreting seismic data in accordance with the 
inventive method described above. 

[0046] While the invention has been described herein with reference to certain examples and 
embodiments, it will be evident that various modifications and changes may be made to the 
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embodiments described above without departing from the scope of the invention as set forth 
in the claims below. 
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